
// separation of variables method

class Orthantst{
	public:
		Orthantst(double *L, double *b, int p, int nu)
		{
			_L=L;
			_b=b;
			_p=p;
			_nu=nu;
		}
		void operator()(int M,double *res)
		{
   			mat u(M,_p);
   			mat y(M,_p);
			double pa,pb,tpz,arg,tp1,ru,at,prod1,sum1,sum2;
   			boost::random::uniform_01<> U;
   			RandomG::Random<boost::random::uniform_01<> > *_runif=new RandomG::Random<boost::random::uniform_01<> >(U);
   			double* prod=new double[M];
			for(int i=0;i<M;i++)
   			{
	  			prod[i]=0;
			}

			for(int j=0;j<_p;j++)
			{
				_st=new boost::math::students_t_distribution<>(_nu+j);
				for(int i=0;i<M;i++)
				{
					if(j!=0)
					{
						sum1=0;
						sum2=_b[j];
						for(int k=0;k<j;k++)
						{
							sum1+=y(i,k)*y(i,k);
							sum2-=_L[k*(_p)+j]*u(i,k);
						}
						at=sqrt((_nu+j)/(_nu+sum1))*sum2/_L[j*(_p)+j];

					}else{
						at=_b[j]/_L[j*(_p)+j];
					}
					pa=cdf(*_st,at);
					pb=1;
    					prod[i]+=log(pb-pa);
					//cout << prod[i] << " ";
					ru = (*_runif)();
		 			arg=ru*pb+(1.-ru)*pa;
					if(arg > .999999999) arg=.999999999;
					if(arg < .0000000001) arg=.0000000001;
					u(i,j) = quantile(*_st,arg);
					prod1=1;
					for(int k=0;k<j;k++) prod1*=(_nu+k+u(i,k)*u(i,k))/(_nu+k+1);
					
					y(i,j)=u(i,j)*sqrt(prod1);
					
				}
				//cout << "Hello"  <<"\n";	
				delete _st;			
			}	
		  *res=prod[0];
		  cout << "res:";
		  for(int k=1;k<M;k++)
		  {
			  if(exp(prod[k])!=0)
			  {
				   *res=log_add(prod[k],*res);
			  	//	cout << *res <<"\n";
			  }
		  }
		   cout << "-logM";
		   *res=*res-log(M);
		   cout << *res <<"\n";
		   delete _runif;
		   delete[] prod;


		}
		mat ui_to_yi(mat ui,int j)
		{
			mat yi=ui;
			for( int i=1;i<j;i++)
			{
				yi(i,0)=ui(i,0)*sqrt(as_scalar(_nu+sum(yi(span(0,i-1),0)))/(_nu+i-1));				
			}			
			return yi;
		}

	private:
		double *_L;
		double *_b;
		int _nu;
		boost::math::students_t_distribution<> *_st;
		int _p;
};
